#############################################################
# 02_run_party_pressure_mexico.R
# This file reads pressure data for Mexican Congress 
#   and estimate Positive-Trait Item Response Model
############################################################
rm(list = ls())

print(version)

# platform       x86_64-apple-darwin17.0     
# arch           x86_64                      
# os             darwin17.0                  
# system         x86_64, darwin17.0          
# status                                     
# major          4                           
# minor          2.1                         
# year           2022                        
# month          06                          
# day            23                          
# svn rev        82513                       
# language       R                           
# version.string R version 4.2.1 (2022-06-23)
# nickname       Funny-Looking Kid           



library (openxlsx)
library (car)
library (data.table)
library (gtools)
library (dplyr)
library (runjags)
library (coda)
library (random)

 
#### Call data #### 
savePathRep <- "/all_data/"

PRI.Pressure<- read.table(file=paste0(savePathRep, "/PRI_pressure_analysis.txt"), sep="\t", header=TRUE)
PAN.Pressure<- read.table(file=paste0(savePathRep, "/PAN_pressure_analysis.txt"), sep="\t", header=TRUE)
PRD.Pressure<- read.table(file=paste0(savePathRep, "/PRD_pressure_analysis.txt"), sep="\t", header=TRUE)

tmp <- list (PRI.Pressure, PAN.Pressure, PRD.Pressure)

 y.count<- array(unlist (tmp), dim=c(nrow(tmp[[1]]), ncol(tmp[[1]]), 3) )



# Number of anchors (for identification)
A <- 2

# Write model
jags.PTIRT <- "model{
## Likelihood
for (i in 1:N){
  for (j in 1:M){
    for (p in 1:P){
      y[i,j,p] ~ dbern(pi[i,j,p])
      probit(pi[i,j,p]) <- beta[j,p]*log(X[i,p]) - log(alpha[j,p])
    }
  }
}
## Priors
for (j in 1:2) {
  for (p in 1:P) {
    beta[j,p] ~ dnorm(0, 0.25)I(0.7,) ##imposing constraints on the discrimination parameters of the first itwo items
  }
}
for (j in 3:(M-A)) {
  for (p in 1:P){
    beta[j,p]  ~ dnorm(0, 0.25)
  }
}
for (j in (M-A+1):M) {
  beta.c[j-(M-A)] ~ dnorm(0, 0.25) ##these items serve as the bridge across parties 

  for (p in 1:P){
    beta[j,p] <- beta.c[j-(M-A)]
  }
}
for (j in 1:M) {
  for (p in 1:P) {
    alpha[j,p] ~ dgamma(.1,.1) 
  }
}
for (i in 1:N) {
  for (p in 1:P) {
    X[i,p] ~ dlnorm(0,10)
  }
}
}"

# Gather data
PTIRT.data <- dump.format(list(y= y.count
                               , N=dim(y.count)[1]
                               , M=dim(y.count)[2]
                               , P=dim(y.count)[3]
                               , A=A
))

# Gather parameters
PTIRT.parameters = c("X", "alpha", "beta", "deviance")

# Initialize parameters
PTIRT.inits.1 <- function() {
  dump.format(
    list(
        X = matrix (data = c(rlnorm(dim(y.count)[1]*dim(y.count)[3], 0, 1))
                             , ncol=dim(y.count)[3]
                             , byrow=T)
      #, X.anchor = rlnorm (1)
      , alpha = matrix (data = rep(rgamma(dim(y.count)[2], .1, .1), 3),
                        ncol = 3)
      , beta = matrix (data = c(rep(rep(1, 3),2), rep(rnorm (ncol(y.count)-A-2), 3), rep(1, 3),rep(-1, 3)),
                       ncol = 3, byrow = T)
      , beta.c=rnorm(A)
      ,.RNG.name="base::Wichmann-Hill"
      ,.RNG.seed=10291
    )
  )
}


# Initialize parameters
PTIRT.inits.2 <- function() {
  dump.format(
    list(
      X =  matrix (data = c(rlnorm(dim(y.count)[1]*dim(y.count)[3], 0, 1))
                   , ncol=dim(y.count)[3]
                   , byrow=T)
      #, X.anchor = rlnorm (1)
      , alpha = matrix (data = rep(rgamma(dim(y.count)[2], .1, .1), 3),
                        ncol = 3)
      , beta = matrix (data = c(rep(rep(1, 3),2), rep(rnorm (ncol(y.count)-A-2), 3), rep(1, 3),rep(-1, 3)),
                       ncol = 3, byrow = T)
      #,beta.c=rnorm(A)
      ,.RNG.name="base::Wichmann-Hill"
      ,.RNG.seed=102011
    )
  )
}

# run model
set.seed(10830)

PTIRT.model <- run.jags(
  model=jags.PTIRT,
  monitor=PTIRT.parameters,
  method="parallel",
  n.chains=2,
  data=PTIRT.data,
  inits=list (PTIRT.inits.1(), PTIRT.inits.2()),
  thin=30, burnin=5000, sample=200,
  check.conv=FALSE, plots=FALSE)

# Gather MCMC chains
chainsPTIRT <- mcmc.list(list (PTIRT.model$mcmc[[1]], PTIRT.model$mcmc[[2]]))

# Check convergence
gD <- gelman.diag (chainsPTIRT, multivariate=F) 
summary (gD[[1]][,1]) 

### Uncomment the next line to save "mexico_ptirt_mcmc.Rda"
# save (chainsPTIRT, file="mexico_ptirt_mcmc.Rda")

